# -*- coding: utf-8 -*-
"""
Created on Sun Apr  2 10:55:34 2017

@author: GangTimes
"""
import matplotlib.mathtext as mathtext
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
import numpy



num=300
def create_data():
    x=numpy.random.uniform(-1,1,num)
    z=list()
    for ix in x:
        temp=numpy.random.random()
        if(temp>0.05):
            if(ix>0.):
                z.append(1)
            else:
                z.append(0)
        else:
            if(ix>0.):
                z.append(0)
            else:
                z.append(1)            
    return (x,numpy.array(z))
def jacobi(theta,x,y):
    a=theta[0,0]
    b=theta[0,1]
    temp=numpy.power(numpy.e,b+a*x)
    s1=-1*numpy.sum(numpy.multiply(x,y))+numpy.sum(numpy.multiply(temp/(1+temp),x))
    s2=-1*numpy.sum(y)+numpy.sum(temp/(1+temp))
    return numpy.matrix(numpy.vstack((s1,s2)))

def hessian(theta,x,y):
    a=theta[0,0]
    b=theta[0,1]
    t1=numpy.power(numpy.e,b+a*x)
    t2=numpy.power(numpy.e,2*b+2*a*x)
    z1=1+t1
    z2=numpy.power(z1,2)
    x2=numpy.power(x,2)    
    s11=numpy.sum((numpy.multiply(t2,x2)/z2)-numpy.multiply(t1,x2)/z1)
    s12=numpy.sum((numpy.multiply(t2,x)/z2)-numpy.multiply(t1,x)/z1)
    s21=numpy.sum((numpy.multiply(t2,x)/z2)-numpy.multiply(t1,x)/z1)
    s22=numpy.sum(t2/z2-t1/z1)
    return -1*numpy.matrix(numpy.vstack((numpy.hstack((s11,s12)),numpy.hstack((s21,s22)))))

def lfun(theta,x,y):
    a=theta[0,0]
    b=theta[0,1]
    return -1*numpy.sum(numpy.multiply(a*x+b,y))+numpy.sum(numpy.log(1+numpy.power(numpy.e,a*x+b)))

def step(theta,delta,x,y):
    alpha=0.2
    beta=0.8
    t=1
    while(True):
        fvall=lfun(theta+t*delta,x,y)
        fvalr=lfun(theta,x,y)+alpha*t*jacobi(theta,x,y).T*delta      
        if(fvall>fvalr[0,0]):
            t=beta*t
        else:
            break
        
    return t

def logistic():
    (x,z)=create_data()
    x=numpy.matrix(x)
    z=numpy.matrix(z)
    error=0.0001
    theta=numpy.matrix([0.5,0.5])

    for i in range(50):
        g=jacobi(theta,x,z)
        h=hessian(theta,x,z)
        
        delta=-(h.I)*g
        decrease=(g.T)*(h.I)*g
        if(abs(decrease[0,0])<error):
            break
        t=step(theta,delta,x,z)
        print(decrease)
        theta=theta+t*delta.T

    print(theta)
    a=theta[0,0]
    b=theta[0,1]
    #x=numpy.linspace(-200,200,2000)
    #绘图
    zh_font=fm.FontProperties(fname=r"c:\windows\fonts\simsun.ttf",size=18)#路径用小写  大写不行 有这个不能保存为eps
    en_font=fm.FontProperties(family='monospace',size=18)
    zf_font=fm.FontProperties(family='Times New Roman',style='italic',size=22)
    fig=plt.figure(figsize=(7,5))
    
    ax=fig.add_axes([0.14,0.14,0.85,0.85])

    plt.rcParams['axes.unicode_minus']=False #用来正常显示负
    ax.set_xlabel('x',fontproperties=zf_font)
    ax.set_ylabel('y or p',fontproperties=zf_font)
    beginpoint=-1
    endpoint=1
    numpoint=5
    ix=numpy.linspace(-1.0,1.0,200)
    iy=numpy.power(numpy.e,a*ix+b)/(1+numpy.power(numpy.e,a*ix+b))
    x= x.A1
    z= z.A1
    line_1=ax.plot(ix,iy,"-r",label=r"y")
    line_2=ax.plot(x,z,"ok",label="Logistic曲线")
    xticklist=numpy.linspace(beginpoint,endpoint,numpoint)
    plt.xticks(fontproperties=en_font)
    plt.yticks(fontproperties=en_font)
    plt.xticks(xticklist)
    ax.legend(loc='best',prop=zh_font)
    fig.savefig('images\logistic2d.pdf',format='pdf')
    plt.show()
    
    

if __name__=='__main__':
    logistic()